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ABSTRACT 

We extend our analysis of the dynamical evolution of simple star cluster models, 
^ . in order to provide comparison standards that will aid in interpreting the results 

' of more complex realistic simulations. We augment our previous primordial-binary 

simulations by introducing a tidal field, and starting with King models of different 
central concentrations. We present the results of iV-body calculations of the evolution 
CM ' of equal-mass models, starting with primordial binary fractions of - 100 %, and N 

, values from 512 to 16384. We also attempt to extrapolate some of our results to the 

larger number of particles that are necessary to model globular clusters. 

We characterize the steady-state 'deuterium main sequence' phase in which pri- 
mordial binaries are depleted in the core in the process of 'gravitationally burning'. 
C r In this phase we find that the ratio of the core to half-mass radius, r c /r^, is similar 

to that measured for isolated systems (Heggie et al. 2005). In addition to the gener- 
ation of energy due to hardening and depletion of the primordial binary population, 
the overall evolution of the star clusters is driven by a competing process: the tidal 
, dissolution of the system. If the primordial binary fraction is greater than 5% and 

the total number of particles TV > 8192, we find that primordial binaries are not fully 
depleted before tidal dissolution, in systems initially described by a King model with 
, a self-consistent tidal field. 

■ We compare our findings, obtained by means of direct iV-body simulations but 

scaled, where possible, to larger AT, with similar studies carried out by means of Monte 
Carlo methods (Fregeau et al. 2003, 2005). We find significant qualitative and quan- 
titative differences with the results in the earlier paper. Some of these differences 
are explicable by the different treatment of the tidal field in the two studies. Others, 
however, confirm the conclusion of Fregeau et al (2005) that the efficiency of binary 
burning in the earlier Monte Carlo runs was too high. There remain unexplained dif- 
ferences, however. In particular, the binary population appears to be depleted too 
quickly, even in the most recent Monte Carlo results. 

Key words: globular clusters - methods: iV-body simulations - stellar dynamics. 
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1 INTRODUCTION 

Several observational surveys of globular clusters have 
highlighted the presence of binaries, which are found 
to consti tute up to 50% of the mass of the core 
fe.g.. see iRubenstein fc Bailvnl Il997t lAlbrow et alJ 1200 it 
iBellazzini et al.l 12001 iPulone et all l2003f) . Such a high 
number of objects is extremely unlikely to have been 
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formed dynamically, as, at the density required to cre- 
ate a significant binary population due to tidal capture 
and three body encounters, the most probable outcome 
from these processes is a merger JChernoff fc Huandll993 
iPorteeies Zwart fc McMillanll2002^ . which could eventually 
lead to the formation of a centra l intermediate mass black 
hole iPorteeies Zwart et al.ll2004l) . Therefore these binaries 
are regarded as being primordial. The primordial binary 
fraction appears to be a key parameter for determining the 
dynamical evolution of a star cluster, as binaries are an ef- 
ficient heating source which can halt core collapse. In addi- 
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tion, interesting astronomical objects can be produced due 
to three- and four-body encounters, such as blue stragglers, 
X-ray binaries and binaries containing pulsars. 

The evolution of star clusters with primordial binaries 
still presents open issues, mainly due to the intrinsic com- 
plexity of numerical simulations of a system where the lo- 
cal dynamical timescale may be many orders of magnitude 
smaller than the global relaxation timescale. (Hard binaries 
have an orbital period of a few hours, while the half mass 
relaxation time t r u may be of order a few billion years.) 
Numerical simulations have thus been performed either us- 
ing approximate algorithms such a s Fokker Planck or Monte 
Carlo methods <Gao et alJ Il99ll : iGiersz fc Spurzeml Eooot 
IFregeau et ail 2003) , which often have to rely on a knowledge 
of the relevant interaction cross sections; or using direct N- 
body codes but emp loying only a modest number of particles 
of order N sa 10 3 (McMillan . et alJll99Ct iMcMillan fc Hud 
Il994l : iHeggie fc Aarsethlll992l) . Both approaches could po- 
tentially lead to misleading results. On the one hand Monte 
Carlo methods rely on input physics that may not accu- 
rately reflect the realistic interactions in a densely populated 
and fluctuating core (e.g. the velocity distribution may be 
assumed to be isotropic), while interaction cross sections in 
the case of unequal masses are not well known. On the other 
hand, extrapolation from direct simulations may prove to be 
non trivial, as it is not clear how to scale the results obtained 
with the number of particles. 

We have thus recently started iHeggie et alJl2005l here- 
after Paper I) a program to survey the basic properties of 
the evolution of an idealised star cluster with a population of 
primordial binaries, with the aim of placing a stepping stone 
between simplified models, such as those of Fokker-Planck 
type, and realistic complex n umerical simulations, such as 
those recently performed by (Porte gies Zwart fc McMillan! 
( 2004). In the first paper of our series we have focused on the 
evolution of isolated, equal-mass models, with a primordial 
binary fraction in the range — 100 % and with a number of 
particles N in the range from 256 to 16384. We have com- 
pared our results not on ly with those of the study carried 
out bv lGao et alJ (Il99 it) (who u sed a Fokker Planck code) , 
and also to some extent with IGiersz fc Spurzeml (I2000T) 
and with iFregeau et alJ (l2003h , but also with a theoreti- 
cal model which pred i cts th e evolution of the core radius 
llVesperini fc Chernofj IT994h . We have shown that signif- 
icant differences arise between direct iV-body s i mulat ions 
and the Fokker Planck calculations of lGao et alJ il99ll) . In 
particular we found that, starting from the same initial con- 
ditions, the number of binaries in the core after the core 
collapse is significa ntly lower in our simulations than the 
number reported bv lGao et al.l (ll99lT) . This could have im- 
portant consequences if one attempted to infer the origi- 
nal primordial binary fraction by observations of the cur- 
rent number of binaries in the cores of globular clusters. 
Another interesting finding is that the ratio of core radius 
r c to half mass radius in the post collapse phase seems 
not to follow the the o retica l expectation from the work of 
IVesperini fc ChernofJ jl994l) . as, by varying the number of 
particles N, we observed a steeper decrease of r c /r^ than 
that expected (oc 1/ log (OAN)). 

In this paper we extend our simulations to include the 
effect of the tidal field of the galaxy. We consider the evo- 
lution of King models with different values of N, the con- 



centration index Wo and the primordial binary ratio /. We 
compare our direct simulations to the Monte Carlo simu- 
lations performed by Fregeau et al (2003, 2005), hereafter 
Fregeau et al. They studied the evolution of isolated and 
tidally truncated star clusters with a population of primor- 
dial binaries in the range 2 — 20 % and a realistic particle 
number (N — 3 10 5 ). Their two-dimensional Monte Carlo 
Method is expected to offer s ignificant advant ages over the 
one-dimensional code used bv lGao et al] (Il99lf) as the influ- 
ence of anisotropy in the velocities can be taken into account, 
and it is unnecessary to assume that the distributions of en- 
ergy and position (for binaries) are independent. The work of 
Fregeau et al. led to some important conclusions. For tidally 
limited clusters they noted for the first time the possibility 
of an initial expansion of the core radius of the cluster in the 
presence of a significant population of primordial binaries, 
when starting from models with a high central concentra- 
tion (i.e. King profiles with Wo > 7). Fregeau et al. also 
showed that, in general, primordial binaries delay the deep 
core collapse phase that is observed in clusters with only sin- 
gle stars, so that the system can be tidally dissolved before 
collapsing. The comparison of our results with their simu- 
lations is mainly focused on fundamental quantities like the 
disruption rate of the binaries, the core radius and the dis- 
solution time-scale. It must be noted, however, that Fregeau 
et al improved their code between the 2003 and 2005 papers 
by the inclusion of direct numerica l integration of encoun- 
ters. They note (IFregeau et al J2005ll that this alters some of 
their earlier results, as we menti on at appropriate poi nts in 
the present paper. In particular IFregeau et al] |2005) note 
that binary burning was too efficient in their earlier models. 

The paper is organized as follows: in the next Section 
we present the parameters of our simulation dataset and the 
numerical method used. In Sec. 3 we give a physical picture 
of the evolution of a star cluster with primordial binaries 
and tidal field that will guide our interpretation of the sim- 
ulations. In Sec. 01and 5 we present our results for runs with 
a tidal field and a tidal cutoff, respectively, including com- 
parison with the results of Fregeau et al. The last section 
of the paper provides a summary and discussions. 



2 SIMULATIONS: SETUP AND ANALYSIS 
2.1 Specification of the models 

The models considered in this paper are tidally limited, with 
stars of equal mass m. The initial distribution is a King 
model with scaled central potential Wo = 3, 7, 11. Our stan- 
dard models have a primordial binary ratio of 10 — 20 %, 
although we have also performed some runs with fewer or 
no primordial binaries (0, 2, 5%) as well as with higher ra- 
tios, 50 % and 100 %. As in Paper I we define the primordial 
binary fraction as: 

f = n b /(n s + n b ) (1) 

with n 3 and n b being the number of singles and binaries 
respectively. This implies that the fraction of the total mass 
in binary stars, in the case of equal masses, is larger in the 
following way: 

f m = 2n b /(n s + 2n b ) 
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For example, for a run with 10% primordial binaries, / = 
10% whereas f m = 2/11 ~ 18%. 

All our results are presented using standard units 
jHeggie fc Mathieul H986'l in which 



G = M 



-4E T 



where G is the gravitational constant, M is the total mass, 
and Et is the total energy of the system of bound objects. 
In other words, Et does not include the internal binding 
energy of the binaries, only the kinetic energy of their centre- 
of-mass motion and the potential energy contribution where 
each binary is considered to be a point mass. We denote the 
corresponding unit of time by 



U = GM 5 ' 2 /(-AE T f /2 



(2) 



in general. For the relaxation time, we use the following 
expression 



3/2 



t r h 



0.138JVr. 
In (0.117V) 



(3) 



where ru is the half-mass radius, and TV denotes the number 
of original objects, i.e. N = n s + rib- When we discuss a 
run with N = 4096 and 50% primordial binaries, we are 
dealing with 6144 stars. We have considered runs with TV in 
the range 512 - 16384. 

Our runs follow the evolution of the system until the 
number of particles drops below 10. However, to avoid possi- 
ble biases introduced by long-living small- TV configurations, 
we define the dissolution time tdis, as the moment when 
98 % of the initial mass of the system is no longer bound. 

Following Fregeau et al., our initial distribution of the 
binaries' binding energies is flat in log scale in the range from 
tmin to 133e m i„, with e min = ma c (0) 2 . Here cr c (0) is the ini- 
tial central velocity dispersion and this choice, if applied to 
an isolated Plummer model, corresponds approximately to 
the standard range adopted in Paper I, i.e. to e m in ~ 5.1 kT 
and 133e m i„ w 680 kT, where the mean kinetic energy per 
particle (over the entire system, with binaries replaced by 
their centres of mass) is 3kT/2. 

To follow the e volution of th e system we have used 
Aarseth's NBODY6 jAarsethl2Q 03f), which has been slightly 
modified to provide additional runtime diagnostics on the 
spatial distributions of single and binary stars. A King model 
with Wo = 7 required almost one month of CPU time on a 
Pentium4 3Ghz PC before tidal dissolution for N = 16384 
and / = 10%; the use of GRAPE hardware would not 
help significantly, as the computational bottleneck is in the 
treatment of the loc al dynamics of binaries (however see 
Maki no"fc Hudli990t for a discussion on the scaling of the 
computational time with the number of particles in the pres- 
ence of primordial binaries). 

The galactic tidal field is treated in our standard sam- 
ple of simulations as that due to a point mass, and the tidal 
force acting on each particle is computed using a linear ap- 
proximation o f the field. The tidal radius r t is defined as 
JSpitzedll987l) : 
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M 



SMgc 



(4) 



where Meal is the galaxy mass, and R,Gai is the distance of 
the centre of the star cluster from the centre of the galaxy. 



We assume Rgo,i is constant, i.e. the cluster describes a cir- 
cular orbit. Thus the equations of motion of a star in the 
cluster are 



x + 2uiy + (k 
y — 2lox 

z 



4u 2 )x = 



+ v z = 



(5) 



where the axes are orientated conventionally with respect 
to the direction to the centre of the galaxy and the galac- 
tic motion of the star cluster, ui is its angular velocity, k is 
its epicyclic frequency, v is its vertical frequency, and the 
right sides are the components of the gravitational accel- 
eration due to the other stars. The value of rt(0) is set at 
the beginning of the simulation to match the cut-off given 
by the adopted initial King model. In our standard runs 
particles are removed from the system when they reach a 
distance from the density centre greater than twice the in- 
stantaneous tidal radius rt, and the value of rt is updated 
during the simulation according to the decrease of M, the 
total mass of the system. 

Note that this method of treating the influence of the 
galactic field is different from the one adopted by Fregeau 
et al., as they imposed a radial "cutoff", i.e. they did not 
consider the tidal force acting on individual particles, and 
simply removed particles wh ose distance excee ded rt (see the 
description of their code in iJoshi et al.ll2000D . At variance 
with our tidal field treatment, the simulations with tidal cut- 
off do not include the Coriolis and centrifugal contributions 
to the particles acceleration. Therefore the two approaches 
have important physical differences. The effects of this dif- 
ferent treatment are discussed in Sec. |5J where we describe 
a subset of simulations (Wo = 7 and Wo = 3, / = 10%, with 
N from 512 to 16384) using the tidal radial cutoff treatment 
adopted by Fregeau et al. 

Our results, unless otherwise noted, are presented by 
applying a triangular smoothing filter with width 2.5 t r h{0) 
(see Fig. 3 and Sec. 3.2 in Paper I for further details). 

A summary of the runs is presented in Tabic 



2.2 Core radius: definition 

In Paper I we adopted the following definition for the core 
radius: 



Ei=l,iV r iPi 



--1.N Pi 



(6) 



where the sum is made over the particles within the half 
mass radius and the density pi around each particle is 
computed from the dist ance to the fifth nearest neighbour 
JCasertano fc Hudll985l) . This definition presents some sys- 
tematic differences (see Fig. from the one adopted by 
Fregeau et al.: 



3cr r 2 



AnGpo 



(7) 



wher e a 2 is the central velocity dispersion (mass weighted 
as in|V csperini fc Chernofllll994|) and po the central density 
JSpitzeJIigslT For a proper comparison between our set of 
simulations and the runs discussed by Fregeau et al. it is 
important to refer to the same quantities, and in this paper 
we adopt as definition for the core radius r c the choice taken 
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Figure 1. Core radius measured using the method in Paper 
I (Eq. El and in Fregeau et al. (Eq. for runs with N = 16384 
and / = 10%. Note that the density-averaged definition of the 
core radius (Paper I) underestimates the core radius at the be- 
ginning of the simulation compared with the theoretical value for 
the corresponding King model (i.e. 0.67 and 0.21 for Wo = 3 and 
7, respectively). In this paper we adopt the definition used by 
Fregeau et al. for the core radius. 

by Fregeau et al., i.e. Eq. Q. We estimated the velocity dis- 
persion o~\ using all the particles within the radius denned 
by Eq. © , while po has been obtained from the Lagrangian 
radius containing 1 % of the instantaneous mass of the sys- 
tem. For the purpose of comparison with Paper I, we also 
report in Tabled] the value of the alternative definition for 
the core radius, r c , at the end of the core collapse. 

2.3 Time-scales 

In the presence of primordial binaries and a tidal field, sev- 
eral dynamical processes operate, each with its own time 
scale. Here we summarize them: 

• Dynamical time (td, cf. eq.(j2j): within a small numeri- 
cal factor, this is the typical timescale for a particle to cross 
the system. It also represents the timescale on which approx- 
imate dynamical equilibrium is established in the case of 
initial conditions starting out of equilibrium. (Note however 
that all our simulations start from approximate dynamical 
equilibrium.) 

• Half-mass relaxation time (t r h, cf. eq.©): this is the 
timescale for energy diffusion in phase space, i.e. for estab- 
lishing thermal equilibrium. It is of order N/ log(0.11A?") td- 

• Total-mass relaxation time (t rt ): this is the timescale for 
relaxation of the cluster as a whole and is defined using the 
tidal radius as characteristic length scale for estimating the 
relaxation time. t T t is proportional to t r h (defined above): 
Ut /Uh oc (r t /r h ) 3/2 . 

• Time for initial core contraction (t cc ): this is the 
timescale on which the core of the cluster contracts (or ex- 



pands) to reach a quasi-steady configuration in which the en- 
ergy released in the core fuels the expansion of the half mass 
radius. The typical value of t cc is of the order of 10t r h(0)- 

• Tidal dissolution time (tdis)'- Time needed to lose 98% 
of the initial mass of the system due to escape across the 
tidal boundary. Depending on the tidal field strength, the 
initial concentration of the system and the initial number of 
particles, this timescale ranges from a few relaxation times 
to a few tens of relaxation times. 

• Tidal half-mass dissolution time (tHaif)- Similar to 
tdis, but defined as the time needed for 50% of the mass 
of the system to escape. 

• Binary burning time (tbb) Time needed to deplete 80% 
of the initial number of binaries. Around this time gravother- 
mal oscillations set in, provided that N is large enough. For 
initial binary fractions above 10%, tbb is of the order of sev- 
eral tens to hundreds of initial relaxation times. 



3 THE PHYSICAL PICTURE: PRIMORDIAL 
BINARIES AND THE TIDAL FIELD 

The evolution of an isolated star cluster made of equal mass 
stars is driven initially by the gravothermal instability. Two- 
body encounters drive a heat flow from the central region 
of the star cluster, which behaves like a self gravitating sys- 
tem with negative specific heat, to the halo. This triggers a 
thermal "collapse" on the timescale of the heat flow process. 

The "collapse" phase lasts until an efficient form of en- 
ergy generation in the centre can stop the process by pro- 
viding an energy production rate equal to the energy loss 
rate by two body relaxation from the region of the core. In 
the case of systems consisting of single stars only, the col- 
lapse is stopped only after one or more binaries are formed 
due to three body encounters in the dense central region. 
The collapse typically takes a time t cc of order 15 t r h(0), 
where t r h(0) is the initial half-mass relaxation time. At this 
point, i.e. at the end of core collapse (sometimes called "core 
bounce"), the central density can exceed 10 4 times the av- 
erage density inside the half-mass radius (for N = 16384, 
or far more for larger N). After core bounce, the generation 
of energy in the core is accompanied by a steady expan- 
sion of the half mass radius. If N is large enough, however, 
binary activity in the core may cause a temporary tem- 
perature inversion which drives the gravothermal instabil- 
ity in reverse. This suppresses binary activity, and further 
collapses may recur, in a succession of "gravothermal oscil- 
lations" (Sugi moto fc Bettwieserl ll983). 

This picture is strongly modified when a population of 
primordial binaries is present, as the generation of energy 
due to existing binaries is more efficient than that due to 
dynamically formed pairs (see Sec. 2 of Paper I). For this 
reason, if the primordial binary fraction is above a few per- 
cent, the initial contraction of the core is more gentle, and it 
is halted at a lower core density (and a larger core radius). 
The timescale for initial core contraction (t cc ) remains of the 
same order of magnitude as in the case where no binaries are 
initially present (see Fig. 17 in Paper I). From then until 
the time when the population of binaries becomes heavily 
depleted, the evolution of the system proceeds with an al- 
most constant ratio of core radius to half mass radius. This 
phase lasts for a time proportional to the initial binary ratio, 



Star Cluster Evolution with Primordial Binaries II 5 



and thereafter a phase of gravothermal oscillations sets in, 
much as for the case where there are no primordial binaries 
(see Figs. 5-6 in Paper I). For an isolated system, provided 
that TV is large enough (i.e. N > 10 4 ), the onset of gravother- 
mal oscillations occurs when the mass of binaries has fallen 
to roughly 10-15% of its initial value flFregeau et aljEooA 
Fig. 4; Paper I, Figs. 5,7). These are empirical results, how- 
ever, and may change for primordial binary fractions smaller 
than the lowest used in these simulations, i.e. 2%. 

The presence of a tidal field introduces a new timescale 
for the evolution of the system, i.e. tdis- In the presence of a 
tidal field, the star cluster steadily loses mass at a rate that 
is almost independent of the properties of the central parts 
of the system, where binaries accumulate (because of mass 
segregation) and act as an energy source for the system. In 
principle the tidal dissolution can be so fast (tdis < tbb) that 
the system does not have enough time to deplete its reser- 
voir of primordial binaries before dissolving. On the other 
hand, if the initial primordial binary fraction is low or if the 
tidal field is weak, then the system can burn almost all of 
the binaries before being dissolved, and it then undergoes a 
phase of gravothermal oscillations, as already described for 
an isolated system, provided that TV is large enough. 

In an isolated system with an initial binary ratio of 10%, 
approximately 100 t r h(0) is the typical timescale tbb required 
for disrupting 80% of the binary population considered in 
our runs (e.g., see Fig. 7 in Paper I). Note, however, that 
the range of binding energy in these runs is 5 — 700 kT, where 
3/2NkT is the total kinetic energy of the cluster (when the 
binaries are replaced by their barycentres) . The time to dis- 
rupt most of the binaries is certainly dependent on the initial 
distribution. 

The time to disrupt most of the binaries can also be 
estimated with the following argument. A primordial binary 
in the range of 5 — 700 kT gives off, on average, around 
100 kT before leavi ng the system (or being dissolved; see 
e.g.. lHut et al Thus if we start with a 10% primordial 

binary population each single star receives around 10 kT. 
The energy loss at the half mass radius due to the tempera- 
ture gradient is of the order of 0.1 kT per star per relaxation 
time, so we have to wait an order of 100 relaxation times be- 
fore depleting the binaries. The same argument suggests that 
the time taken to burn most of the binaries is proportional 
to /. 



4 TIDAL FIELD RUNS: RESULTS AND 
COMPARISONS 

In this Section we discuss the properties of our simulations, 
our main purpose being to set out the essential empirical 
facts about the evolution of systems with primordial bina- 
ries in a galactic tidal field. We shall also compare our results 
with those obtained by Fregeau et al., though they treated 
the tidal field as a cutoff. To elucidate this comparison, we 
were motivated to carry out further JV-body runs with the 
same treatment (i.e. a cutoff). We defer presentation of those 
calculations to the next section, however; in the present sec- 
tion we consider only runs that include our most realistic 
treatment of the tidal force (eq.@). 



4.1 Total mass and dissolution timescale 

Our results are summarized by the values of tdis in Tab. [3] 
and Tab|21 and details for a subset of the runs are presented 
in the top panel in Figs. 12151 these are intended also to facil- 
itate comparison with corresponding plots in Fregeau et al., 
as specified in the captions. The evolution of the total mass 
in our different runs depends on three initial parameters: (i) 
the initial fraction of primordial binaries, /; (ii) the scaled 
central potential of the initial King model, Wo; and (iii) the 
initial number of objects, N. We consider the dependence 
on each in turn. 

4-1.1 Dependence on f 

The initial fraction of primordial binaries has only a small 
influence on the rate of mass loss of the system and on tdis- 
The dependence is, however, systematic, as can be seen from 
Fig. [(J (Though this figure shows runs with Wo = 7, the 
trend shown is representative of all concentration indexes 
that we have studied.) The simulation with the longest dis- 
solution time is the one starting with 100% of primordial 
binaries. This may not be surprising, as this run has the 
highest total number of stars and so, after disruption of a 
given proportion of binaries, it has the longest relaxation 
time. On the other hand it is closely followed by the simu- 
lations with low binary content (2 % and 5 %) and then by 
the run with single stars only, which is dissolved in 90% of 
the time needed in the 100 % case. This may be due to the 
fact that simulations with / < 5% deplete their primordial 
binaries before dissolution (tbb < tdis); this leads to a deep 
core collapse, which creates a strongly bound core which is 
rather resilient to dissolution. Curiously, simulations with 
20% of primordial binaries are subject to the fastest tidal 
dissolution, and this could be due to effects related to mass 
segregation, which is absent if / = or 1, and which tends 
to expel single stars. 

4-1.2 Dependence on Wo 

Consider, as an example, runs with N = 4096 and / = 10%. 
A low concentration King model with Wo = 3 is tidally dis- 
solved in approximately 13 t r h(0) (Table 01 , while the sys- 
tem can survive for approximately 50 t r h(0) for Wo ~ 7; but 
Wo = 11 is intermediate between the two (tdis ~ 25 t r ft(0)). 
This non-monotonic dependence on concentration index Wo 
is a general feature, independent of the number of particles. 
It may be understood in terms of the non-monotonic depen- 
dence of the tidal radius on Wo, at fixed total mass an d total 
energy (Tab[T] see also Fig. 8.3 in lHeggie fc Hut 120031) . (An- 
other way of looking at this is that, in the units we adopt, 
the strength of the tidal field varies non-monotonically with 
Wo.) Though the difference in rt between the cases Wo = 7 
and Wo = 11 may not seem great, two further factors may 
be relevant: (i) rn also increases slightly between Wo = 7 
and Wo = 11 (Table [TJ, which would increase the vulnera- 
bility of the more concentrated system, even if the tidal field 
were the same; (ii) high-concentration models begin with a 
short period of core expansion (Fig|^J, and the associated 
burst of energy generation enhances tidal overflow. 

If the tidal dissolution time is expressed in units of the 
total-mass relaxation time (t r t), then the differences be- 
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Table 1. Initial tidal and half-mass radius rt(0), r'h(O) (in natural 
units) as functions of the scaled central potential Wo 

Wo r t (0) r h (0) 

3 3.15 0.84 
7 7.00 0.81 
11 6.62 1.01 



tween the different King models are greatly reduced. By 
accounting for the factor (rt/rh) 3//2 , so that the time is in 
units of trt, the difference in tidal dissolution time between 
the various King models considered here is within 20% for 
TV = 4096 and within 10% for TV = 16384. 



4-1.3 Dependence on N 

In our set of simulations we observe a marked dependence 
of the dissolution time-scale tau ° n the number of particles 
used (see Tab. [3] and Fig. 0. This is not surprising, as a 
similar effect is well known in simulations with the tidal 
force modeled as in our eq. though using single particles 
only: the time for the loss of half of the mass of a cluster is 
approximately proportional to i r h(0) 3 ^ 4 , i.e. 

/ In (0.117V) \ 1/4 , . 

tBalf OC " J *rfc(0), (8) 

(Baumgardt 2001). Besides empirical support from TV-body 
simulations, this formula can be understood on the basis 
of the combined effect of the tidal field and of two body 
relaxation (Baumgardt 2001). We find (FigEJ that it also 
provides a reasonable fit to the dissolution time tdis of mod- 
els with primordial binaries (for the parameters specified in 
the caption of the figure), though the observed dependence 
on TV is slightly stronger. Therefore we can take advantage 
of eq.lJSJ for the approximate extrapolation of results on the 
dissolution time tdis to larger num bers of parti c les TV . Note, 
however, that eq.(|HJ is derived in iBaumgardtJ i200ll) using 
simulations with TV ^ 16384; it is expected to be valid up 
to TV < 10 6 at most, as for larger TV the equation eventually 
leads to tualf < t r h(0), a clearly unphysical result. 

Assuming, then, that this scaling applies to the disso- 
lution time tdis, we can compare our results with those in 
Fregeau et al. If we extrapolate to TV = 3 10 5 , we find that 
tdis decreases by a factor « 2.6 from the results obtained for 
the simulations with TV = 4096 and by a factor m 1.9 from 
the runs with TV = 16384. This would mean that a simula- 
tion starting from Wo = 7 with 10 % binaries and TV = 3 10 5 
would take w 18 t r h(0) to dissolve. For comparison Fregeau 
et al. measure tdis > 38t r h(0) 1 for these initial conditions. 
This and other comparisons are given in Table H (second 
and fourth columns). An understanding of such differences 
requires consideration of the treatment of the tide, and so 
we defer further discussion to Sec^l 
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Figure 2. Time dependence of the total mass and of the mass 
in binaries (upper panel) and of the half mass and core radius 
(bottom panel). The dotted line is the half mass radius for singles 
while the dashed line in the lower panel is the half mass radius 
for binaries. The lowest curve is the core radius of the system (in 
units of the initial half mass radius). The simulation has been 
performed with 16384 particles and 10% of binaries starting fr om 
W = 3. It is the equivalent of Fig 13 in lFregeau et al] <2003l) 
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Figure 3. Like Fig.[5]for a simulation starting from Wo = 7 with 
16384 parti cles and 10% of binaries. It is the equivalent of Fig 10 
in lFregeau et al] ]2003l) 



1 Note that our definition of t r h(0) differs by a factor 
ln(0.11TV)/ln(0.4iV) from that used by Fregeau et al.. This has 
been taken into account in this paper to ensure a proper compar- 
ison. 
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Figure 4. Like Fig.|2]for a simulation starting from Wo = 7 with 
16384 parti cles and 20% of binaries. I t is the equival e nt of Fig 11 
in IFregeau et alJ i2003fi or Fig.3 in IFregeau et al.l J2005D . The 
results from a new, unpublished run by Fregeau (Fregeau 2006, 
private communication) are overplotted in red. 
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Figure 5. Like Fig. |2] for a simulation starting from Wo = 11 
with 163 84 particles and 10% of binaries. It is the equivalent of 
Fig 12 in IFregeau et alJ <2003h 



4.2 Core and half mass radius 

In the presence of a tidal field, the system evolves toward 
conditions in the centre very similar to those that we ob- 
served for isolated systems. The core radius evolves, usually 
in a few t r h, to reach a ratio r c /rh that is close to that 
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Figure 6. Evolution of the total mass for a King models with 
Wq = 7 and different ratios of primordial binaries. The simula- 
tions have been performed with 4096 particles. 



Table 2. Comparative results on the dissolution timescale 
tdis/t i,(0) f° r simulations with / = 10 % 

Wo THH field THH cutoff Fregeau et al. cutoff 

3 5 15 13 

7 18 20-40 >38 

11 11 N/A 27 



Notes: The values for t^is for our runs (THH) have been extrap- 
olated to N = 3 ■ 10 from N = 16384 (see discussion in the 
text). In the second column we give t^ ia for our standard tidal 
field runs, extrapolated by using eq. l8l . In the third column we 
report our results from runs with a tidal cut-off, with constant 
extrapolation for Wo = 3 and a range of values for Wo = 7 (see 
Sec|FJ. The fourth column contains the results from Fregeau et 
al. 

attained at the end of core collapse during the evolution 
of isolated Plummer models (see Sec. 4.1 of Paper I) with 
a similar fraction of primordial binaries. As is depicted in 
Fig. |H1 the evolution of the core to half mass radius, after a 
first adjustment phase, is largely independent of the initial 
conditions considered. Thereafter this ratio is almost con- 
stant, as long as there are binaries to burn or until tidal 
dissolution has reduced the total mass below 10% of the 
initial value. 

We discuss the value of the core radius in further de- 
tail below, but summarise here the behaviour of the half 
mass radius. This depends on the strength of the tidal field, 
which depends on Wo (Sec l4.1|l . For Wo = 7, ru remains 
almost constant, up to the final stages of the life of the sys- 
tem, when it eventually falls to zero. For stronger tidal fields 
(models with Wo ~ 3 and 11) the half mass radius decreases 
steadily. Our interpretation is that, for Wo = 7, the tendency 
toward expansion (due to the energy generated in the core) 
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Figure 7. Tidal dissolution t^is time for models with Wo = 7, 
/ = 10% and a galactic tidal field, for different values of N, 
compared to the result by Fregeau et al. (square point; obtained 
using a tidal cut-off). The points associated with runs using a 
tidal field are compared with t he theoretical scaling (red solid 
line) given by iBaumgardd JsiOOl) for tidal field runs. 



is balanced by the mass loss due to the tidal field, while the 
tide dominates in the other two cases. 

Now we return to the evolution of the core radius, 
and consider first the initial phase of contraction or expan- 
sion. The evolution in the first few relaxation times depends 
on how the initial value of r c compares with the quasi- 
equilibrium value attained in the intermediate steady binary 
burning phase. Thus runs starting initially with a relatively 
big core, such as Wo = 3 (r c (0) w 0.6r^(0), Fig. |2J have a 
rather deep collapse; after about 8t r h(0) in our simulation 
with N = 16384 and / = 10% we measure r c /r h (0) ~ 0.05. 
For Wo = 7, / = 10% and N = 16384 (see Fig. EJ the 
duration of the initial contraction is shorter, and thereafter 
the core radius is somewhat larger. On the other hand runs 
that start with a small, concentrated core (Wo = 11), have 
a first rapid expansion of the core radius (see also Fig. [SJ up 
to r c /r h (0) « 0.09 for our N = 16384 run with / = 10%. 
This can be understood in terms of the energy balance of 
the core: if the density is too high, the binary heating ex- 
ceeds the cooling by heat transport to the halo, so that the 
core has to expand to lower density to reduce the rate of 
energy generation. Incidentally, though we have focused on 
numerical values for runs with / = 10%, we find that the ef- 
fect of primordial binaries saturates by around / = 10% (see 
Fig. 0, much as we have noted in the evolution of isolated 
clusters in Paper I (see also Vesperini & Chernoff 1994). 

Now we turn to the value of the core radius during the 
phase of steady binary burning, and, in particular, its de- 
pendence on the number of particles. Our results are sum- 
marised in Fig. 1 10 1 At va riance with Paper I, here the 
IVesperini fc Chernofll ill 9941) model provides an excellent fit 
to the dependence on iV of the ratio of the core to half mass 



Figure 8. Evolution of the ratio of the half mass to core radius 
(upper panel) and of the concentration parameter c = log rt/r c 
for different King models with / = 10%, 20%. The solid line refers 
to simulations starting from Wo = 7, the dotted line to Wo = 11 
and the dashed line to Wo = 3; the number of particles used 
is 16384. The evolution quickly erases differences due to initial 
conditions, and during binary burning c and r^/r c evolve very 
similarly. This figure can be directly compared with Fig. 16 in 
Fregeau et al. Note that we kept the same range in the y-axis to 
highlight the differences between our runs and theirs. 

radius, provided that in the Coulomb logarithm the value 
0.11 is used for t he coefficient of N. Thus th e formula used 
for the fit is (see lVesperini fc Chernofj|l994l and Paper I): 

_Tc _ a <fa(l — 4>b)^bs + <fbVbb 

r k ~ Iog 10 (0.11JV) (l + <« 4 ' U 

where the various quantities are defined as in Paper I. Thus 

the quantity a is a parameter depending on 7 and v c /vh, 

where 7 is the ratio of the expansion timescale r^/fj, to 

trh, and v c and vu are the velocity dispersion in the core 

and at the half-mass radius, respectively; we have assumed 

typi cal values for these parameters: 7 = 10.5 and Vc /vh = 

y/2 jHeerie fc Aarsetblll992l; iGoodman fc HutJll989l) . The 

quantity <f>b is the binary fraction in the core defined as 

, rib 

<Pb = : • 

n s + rib 

Finally, fii, s and fibb are coefficients for the efficiency of 
binary-single and binary-binary burning and depend on the 
distribution of binding energy of the binaries. They have 
been computed as in Paper I assuming a flat distribution in 
the logarithm of the binding energy from 12 kT to 300 kT. 
a, jj,b 3 and fibb do not significantly depend on the number of 
particles N (see Sec. 5.1.2 in Paper I). 

With regard to the core radius there are many quan- 
titative and even some qualitati ve differences b etwee n our 
results and the earlier results oflFregeau et alJ ll2003l) . But 
they point out iFregeau et al.ll2005t) that the core radius is 
one parameter which has turned out to behave rather differ- 
ently following the code improvements summarised in their 
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Figure 9. Dependence of the ratio r c /r^ on time (in units of 
tr-ft(O)) f° r simulations starting from a King model with Wo = 7 
and different ratios of primordial binaries. The simulations have 
been performed with 4096 particles. 
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Figure 10. Dependence of the ratio r c /r^, averaged for 10 t r h(0) 
after core collapse, on the number of particles. Simulations start 
from a King model with Wo = 7 and / = 10% using a tidal 
field (crosses) or isolated runs (squares). Points associated with 
TV from 512 to 2048 have been obtained by averaging over a 
sample of 16 simulations f or each N. A comparison with the 
IVesperini &; Chernofll ll994l) model is drawn: the agreement (us- 
ing £og(0.11TV) instead of log(0AN)) is excellent for the tidal field 
runs. The adopted values for the parameters used are given in the 
text. 
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Figure 11. Like Fig. |2] for a simulation starting from Wo = 7 
with 4096 particles and 5% of binaries. 

later paper. Indeed the one case for which a direct com- 
parison with their most recent results is possible is shown 
in Fig|l(Wo = 7, / = 20%). Both models (our TV-body 
model and their Monte Carlo model) show an initial con- 
traction followed by a phase of stead y binary burning. The 
comparison with the new results by iFreeeau et al.l <l2005l) 
is quantitatively satisfactory, except that, in the phase of 
steady binary burning, our core radius is about twice as 
large as theirs. However, it is known that the core radius in 
this phase is TV-dependent (see above, and Paper I) , and the 
sense of the difference is consistent with this. 

As already mentioned in Sec El the late evolution of the 
core radius depends on whether the binary population can 
be heavily depleted by the time of dissolution, tdis- The evo- 
lution of the binary population is discussed in the following 
section, and so we defer further discussion until then. 

4.3 Evolution of the primordial binary population 

The destruction of the primordial binary population is de- 
picted in the top panel of Figs. 12151 The number of binaries 
decreases due to both destruction of the softer pairs in the 
centre of the system and the ejection of stars across the tidal 
boundary. The changing relative contribution of these two 
processes can be inferred from Fig. 1121 upper panel. At first 
the ratio of binaries to singles decreases, as at the beginning 
of the simulation relatively many soft binaries (i.e. those 
with Ei, < 15kT) are present. However at later times the 
relative number of binaries to singles starts to rise, as most 
of the surviving binaries are hard to disrupt in three- or four- 
body encounters. In addition tidal ejection is more probable 
for singles, as binaries have sunk toward the core of the clus- 
ter by mass segregation (see bottom panel in Figs. 121 51 and 
Fig.im 

There are some quantitative differences between our re- 
sults and those of Fregeau et al., even after the recent im- 
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provements to their technique. For example they point out 
ilFreeeau et"ai] 120051) that, for a model with Wo = 7 and 
/ = 20%, the binary fraction shrinks to about 10% before 
recovering to about its initial value shortly before tdis- This 
differs considerably from our result in Fig |12l (upper panel). 
Another example is illustrated in the upper panel of Fig^] 
where it is clear that the disruption of binaries in the most 
recent available run by Fregeau et al. appears to be rather 
faster than in our run starting from similar initial conditions. 

Now we consider models in which the population of bi- 
naries can be strongly depleted before dissolution of the clus- 
ter, i.e. ttb < tdis- We have already remarked fSec l4.H that 
tdis increases as TV decreases, and is longest for Wo = 7 (out 
of the three values we adopted). It has also been argued 
(Sec[3J that the time taken to burn most binaries is roughly 
proportional to /. Based on these considerations, it may be 
expected that the only runs that may burn most of their 
binaries within the lifetime of the system are those in which 
/ or TV is sufficiently small, and that intermediate values of 
Wo are most favourable. Indeed, we found that the only runs 
which did so were those with Wo = 7 and TV < 4096 and 
/ < 5%. (Note, however, that Wo = 7 is the only value of 
Wo for which we carried out runs with TV < 4096, see Table 
(21) At the point where the binaries become heavily depleted, 
there is a deep core contraction, similar to what happens in 
runs with single stars only; however, here TV is too low for 
gravothermal oscillations to occur. In Fie llll we depict the 
evolution of a Wo = 7 model starting from / = 5% and 
TV = 4096: by about 30 t r h{0) the fraction of binaries has 
fallen below 20 % of its initial value and the core radius 
contracts sharply. 

Next we comment in particular on the evolution of the 
binary population in the core. At the end of our simulations 
starting from / = 10 — 20% we find that the fraction of bina- 
ries in the core has increased by at least a factor 4. In general 
the density of binaries in the core is in broad agreement with 
the numbers observed in Paper I and does not depend sig- 
nificantly on the details of the initial conditions if / is in 
this range (see bottom panel of Fig. 1121 where runs starting 
from Wo = 7 and Wo = 11 show a rather similar behaviour). 
However there are differences in detail in the time evolution. 
In an isolated cluster the proportion of binaries in the core, 
after the sharp initial increase caused by mass segregation, 
declines slowly; in these tidally limited models, on the other 
hand, there is a slow increase. The Wo = 3 King model is 
an extreme example, as there is an enhanced rate of ejection 
of singles relative to the rate of disruption and ejection of 
binaries, because of the relatively small tidal radius; in this 
situation the binary fraction in the core increases rapidly 
throughout the evolution. 

The evolution of the internal binding energy (see 
Fig. 1141 is similar to that observed for isolated systems, with 
a preferential decrease in the number of less bound binaries. 
Furthermore, as expected, there is a significant correlation 
of the internal binding energy with the radius (see Fig. 1151 . 
After approximately 15 t r h{0) the survival probability in the 
core for binaries with binding energy below « 20 kT is low, 
as they have either been destroyed or hardened. We thus 
see that softer binaries are mainly present around the half 
mass radius and in the halo, while the core is dominated 
by hard binaries (with occasionally some short-lived low en- 
ergy binaries that have been formed dynamically). At later 
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Figure 12. Dependence of the binary fraction on time (in units 
of t r 7,(0)) for simulations starting from King models with Wo = 3 
(dotted lines), Wo = 7 (solid lines), and Wo = 11 (dashed lines) 
for / = 10%, 20%. The upper panel refers to the total number 
of binaries, the lower to the binaries in the core. The simulations 
have been performed with 16384 particles. 

stages of the evolution (see bottom right panel of Fig. I15|l . 
the softer binaries have almost completely disappeared. 

This picture is in qualitative agreeme nt with the results 
already known in the literatur e (see e.g., iMcMillan fc Hull 
Il994l : lOiersz fc Spurzemll200o| and Fregeau et al.). As our 
resolution is limited by the modest number of binaries that 
are present in our runs, it is hard to attempt any quantita- 
tive comparison. 



5 TIDAL CUTOFF RUNS: RESULTS AND 
COMPARISON 

It was pointed out in Sec 14.1.31 that our result on the life- 
time differed significantly (at least in the one case consid- 
ered) from that of Fregeau et al. This conclusion, however, 
required a knowledge of the TV-dependence of the lifetime, 
and it is known that this in turn depends on the treatment 
of the tide, which differs between the work of Fregeau et 
al. and the TV-body simulations which we have presented 
so far. For the case of a tidal cutoff, as a dopted by Fregeau 
et al., it was found bv iBaumgardd <l200ll) that the half-mass 
dissolution time-scale tHalf (expressed in relaxation times) 
does not depend on the number of particles of the system (as 
we assumed for the purpose of our com parison) if TV > 4096. 
Note, however, that iBaumgardt) i200ll) used a fixed tidal ra- 
dius and, as we have mentioned, measured the timescale of 
mass-loss by the half-mass time tHalf, whereas our discus- 
sion has focused on the time for loss of 98% of the mass, 

tdis • 

Because of these complications, for a better comparison 
with Fregeau et al.we have run a second series of simulations, 
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Figure 13. Evolution of the distribution of binary radius in units 
of the instantaneous core radius. The ordinate is the number of 
binaries per bin, which are of length 0.5 in log e . Curves have 
unit normalization, and are shown for t/t r h(0) = 0, 8, 16, 24. The 
dashed line is the position of the half mass radius. The simula- 
tion has been performed with 16384 particles and 20% of binaries 
starting from a King model with Wo = 7. 



Figure 15. Evolution of the distribution of binaries in radius 
and binding energy. The radius is given in units of the instan- 
taneous core radius; the energy in units of kT. The panels refer 
to t/t r h(0) = 0,8,16,24. The dashed line is the position of the 
half mass radius. The simulation has been performed with 16384 
particles and 20% of binaries starting fr om a King model wi th 
Wo = 7. It can be compared to Fig. 15 in Frcg eau et alJ (2003) 
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Figure 14. Evolution of the distribution of binary binding energy 
(given in units of kT). The ordinate is the number of binaries per 
bin, which are of unit length in log 2 . Curves have unit normaliza- 
tion, and are shown for t/t r f l {0) = 0, 8, 16, 24. The simulation has 
been performed with 16384 particles and 20% of binaries starting 
from a King model with Wo = 7. 



in which the tidal field we have used hitherto has been re- 
placed by a tidal cutoff. We have ensured that the cutoff 
radius decreases in the correct manner (eq.Q) as the mass 
inside this radius decreases, exactly as adopted by Fregeau 
et al.The initial parameters of these models are Wo = 7 and 
Wo = 3, / = 10%, with TV in the range 512 to 16384. 

For the series starting from Wo = 3 with / = 10 % 
(i.e. with a strong tidal field, Sec 14. 1.2^ . we do not observe a 
significant TV-dependence of the half-mass dissolution time 
tHalf in units of the relaxation time, except for a modest 
(« 10%) variation from TV = 4096 to TV = 8192. This is 
consiste nt with the experim ents (without primordial bina- 
ries) by iBaumgardtl ll200ll) , despite his use of a constant 
tidal field. Moreover, it is also consistent with unpublished 
results by one of us using the initi al conditions o f the collab- 
orative experiment summarised in iHeggiel 1120031) . except for 
replacement of the tidal field by a tidal cutoff; in this case the 
cutoff was evolved in the correct manner as mass was lost. 
For the Wo = 3 runs with tidal cut-off the agr eement with 
the dissolution time tdis of lFreeeau et al.l ll2003l) is good: we 
extrapolate tdis ~ 15 t r h{0) and they have tdis = 13 t r h(0). 

By contrast, for the series with Wo = 7 and f — 10 % 
and a cutoff, we observe an TV-dependence similar to that 
found when a tidal field is used, as in our standard runs (see 
Fig. 1161 . On the other hand it is not known why the TV- 
dependence should differ for the cases Wo = 3 and 7 when a 
cutoff is used. One possibility is that it is a "small- TV" effect 
which is more important in a system with a small core; if so, 
the TV-dependence of the dissolution time tdis will eventually 
flatten out for larger TV. The presence of primordial binaries 
does not seem to be the issue; we have run a few simulations 
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with a tidal cut-off, starting from a Wo = 7 model and no 
binaries, and find that the half-mass dissolution time tnaif 
is almost equal to that measu red when binaries are present. 
Tani kawa fc Fukushigd i2005t) have also observed differences 
(in the scaling of the half-mass time tuai / with TV) between 
models with Wo = 3 and 7, though their models adopted a 
tidal field, and the potential was softened. They tentatively 
attribute the difference to the time taken for a particle to 
escape. 

Whatever the explanation, this (unexpected) TV depen- 
dence of tdis in our tidal cutoff runs with Wo = 7 makes the 
extrapolation of our results obtained with the tidal cut-off 
prescription highly uncertain. While the dissolution time of 
our model with W = 7, / = 10% and TV = 16384 in a tidal 
cutoff is tdis — 41.2i r k(0) (Table|3Jl, the inferred dissolution 
time for a Wo = 7 / = 10% model with TV = 3 ■ 10 5 may be 
in the range of 20 — 40t r h(0) (Tabled col. 3), depending on 
whether tdis is proportional to t r h (0) or scales in accordance 
with eq.lJSJ. If the actual value is close to the upper limit, 
then the tidal dissolution time tdis in our Wo = 7 cutoff 
models would compare marginally well with Fregeau et al., 
as they find tdis > 38. Furthermore, recent improvements to 
their code appear to have reduced the lifetime slightly for a 
model with Wo = 7, / = 20%, and a similar reduction for 
/ = 10% would improve the agreement between our TV-body 
and their Monte Carlo results. 

We close this section with a brief comparison between 
other properties of our runs with a tidal field and those 
with a tidal cutoff. Because of the limited scope of our runs 
with a tidal cutoff, we are concerned less here with general 
trends; but these results will be useful for comparison with 
simplified treatments in which it is necessary to use a cutoff, 
such as Fokker-Planck models; they also illustrate the kinds 
of results which may be sensitive to the manner in which 
the tidal boundary conditions are treated. 

For our run with Wo = 3, / = 10%, TV = 16384 and 
a cutoff fFig. 1171 . the initial contraction of the core radius 
is slower than that observed in the run with a tidal field 
(see Fig.|5J: after 8t r h(0), r c is approximately twice as large 
in the cutoff run. However just before tidal dissolution the 
values of the core radius in the two runs are quite similar. 
Actually, if the time is normalized by the dissolution time 
tdis, the evolution of r c looks very similar, except at the very 
last moment, where r c decreases for the run with a cutoff. 

By contrast, the run with Wo = 7, / = 10 %, TV = 16384 
and a cutoff (Fig. I18H is much more similar to the corre- 
sponding run with a tidal field (Fig.|3J. The dissolution time 
tdis scales for the two tidal treatments differ by only w 25%, 
and the core radius evolves in a similar way. However the 
half mass radius for the run with a cutoff initially increases 
in much the same way as in the run reported by Fregeau 
et al., while in our run with a tidal field it remains almost 
constant. 



6 DISCUSSION AND CONCLUSIONS 

In this paper we have continued our investigation of the 
dynamical evolution of stellar systems with primordial bi- 
naries. In Paper I we considered simple isolated models 
with equal mass stars and a primordial binary population 
in the range — 100 %. Then the main parameters that 
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Figure 16. Tidal dissolution time tdis f° r a Wo = 7 / = 10% 
model with a tidal cutoff, with different values of N, compared to 
the result by Fregeau et al. (square point; obtained using a tidal 
cut-off). The lower series (3-sided stars) gives the half mass disso- 
lution time tjj a if, compared to some runs starting from Wo = 7 
and no binaries (triangles). The scaling with N of both tdis and 
tHalf appears similar. 
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Figure 17. Like Fig. [5] but for a tidal cutoff simulation starting 
from Wo = 3 with 163 84 particles and 10% of binaries. It is the 
equivalent of Fig 13 in Freg eau et alj $2003 ) 
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Figure 18. Like Fig.EJfor a tidal cutoff simulation starting from 
Wo = 7 with 163 84 particles and 10% of binaries. It is the equiv- 
alent of Fig 10 in lFreeeau et alj fc003ft 

control the evolution of the system are the primordial bi- 
nary fraction / and the number of bodies N. However, for 
/ > 10% the efficiency of energy generation due to three- 
and four-body encounters saturates, so that little difference 
was observed in the size of the core radius for runs starting 
with the same number of particles and / > 10%. This be- 
h aviour had already been pred icted by the theoretical model 
of lVesperini fc Chernofll dl994T> . and our runs confirmed this 
expectation. Interestingly, however, the scaling of the core 
radius with the num ber of particles was not quit e the same 
as that predicted bv lVesperini fc Chernofll J1994T) . 

Here we have added an important new ingredient to our 
simulations: a tidal field. The presence of a tidal field intro- 
duces a new time-scale in the simulation, the time needed for 
the dissolution of the system, which is usually shorter than 
the time needed to deplete the primordial binary population 
for simulations starting from / > 10% and N > 4096. We 
have considered runs starting from King models with initial 
central potential Wo = 3, 7, 11 (and thus a different inten- 
sity of the tidal field, in the units adopted) with to 100 % 
primordial binaries, and N from 512 to 16384. 

By comparison with the results of Paper I, our Fig. |H| 
clearly confirms that the evolution of the core to half mass 
radius proceeds much as for isolated models, after an initial 
transient; the models starting with Wo = 7 and 11 evolve 
very similarly, for initial binary fractions of / = 10%, 20%. 
The King model with Wo = 3, however, nicely illustrates the 
effects of a strong tidal field: the evolution of r c /ru tries to 
go toward the steady value reached by the other runs with 
higher Wo, but the tidal field manages to dissolve the system 
just when this common r c /rh value is about to be attained. 

One objective of this paper was also to compare the 
results obtained by means of our direct iV-body simula- 
tions with the outcome of simulations able to employ re- 
alistic number of particles by using approximate methods, 



like Monte-Carlo and/or Fokker Planck approaches. Just as 
in Paper I we comp ared our work with the milestone study 
of iGao et alJ (ll99lT) . here we took as comparison the recent 
work carried out by Fregeau et al. There are considerabl e 
differences with the results reported in lFreeeau et al] (2003), 
but they have since reported that several results change sig- 
nificantl y following an impr oved treatment of binary inter- 
actions (iFregeau et al.ll2005f) . As a result of these improve- 
ments they have discovered that the emission of energy in 
encounters with binaries was too high in the earlier models. 
Therefore we have focused our comparison on one case which 
they illustrated in some detail in their later paper. While it is 
promising to report that most differences have been cleared 
up by the improvements, there are exceptions. One is that 
the binaries still appear to be depl eted too rapid l y, eve n 
in the most recent runs (FigEl cf. IFregeau et ail i2005ft . 
Fig. 3). There remain uncertainties also about the lifetime, 
but this depends on perplexing problems in the scaling with 
N: we find e viden ce that the tidal cutoff scaling derived by 
iBaumgardtl ]200l]) for single stars with a fixed cutoff radius 
does not seem to be applicable to the self-consistent pre- 
scription of a cutoff radius prescription used by us (for some 
simulations) and by Fregeau et al. Therefore the best ap- 
proach to understanding the remaining differences between 
our runs and those of Fregeau et al. would be to directly 
compare a set of Monte Carlo simulations starting from the 
same initial conditions as we have adopted. 

We finally n ote that, in the presence of a tidal field, the 
prediction of the lVesperini fc Chernofll il994h model, on the 
size of the core as a function of the number of particles used 
in the simulation, is in detailed quantitative agreement with 
the values measured in our set of simulations with Wo = 7. 
This fact comes as a surprise, becau se, as we discussed in 
Paper I, the N dependence given bv lVesperini fc Chernofll 
( 1994) was not verified for isolated clusters: however, the 
model does not make any assumption about the presence or 
absence of a tidal field. The presence of a tidal field seems 
thus to introduce a difference in the scaling with N of the 
ratio rc/rh in the steady burning phase, which it would be 
interesting to understand better. 
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Table 3. A^-body models with Pregeau et al.-like initial conditions. In this table we report summary information for all runs with 
N 4096. For runs with smaller N we report averaged values over 16 simulations for each entry line. The entries at the bottom of the 
table refer to simulations with tidal cut-off as in Fregeau et al. Column entries are: the number of bodies N, the initial value of Wo, the 
primordial binary ratio /, time of core collapse t cc , time of dissolution t dis (computed at the point where only 2% of the initial mass of 
the system is left), core radius at core collapse r c , and core to half mass radius after core collapse r c /r^(t cc ) (averaged for 5 t r h(0) after 
core collapse), fraction of binaries to singles in the core after core collapse (N b /N s ) c (t cc ) (same average as for core to half mass radius), 
and relative mass of the system at core collapse M(t cc )/M(0). 
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Wo 


/ 


tcc 


tdis (98%) 


f 'c (tcc ) 




r c /rh(tcc) 


(N b /N s )c(tcc) 


M(t cc )/M(0) 


512 


7.0 


10 


1.6 


81.2 


0.147 


0.161 


0.179 


0.35 


0.99 (averaged) 


1024 


7.0 


10 


4.0 


66.8 


0.129 


0.147 


0.157 


0.44 


0.97 (averaged) 
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10 
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0.143 
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10 
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0.15 


4096 


3.0 


10 


10.5 


13.0 


0.073 


0.080 


0.175 


2.53 


0.17 


4096 


3.0 


20 


9.5 


13.4 


0.085 


0.097 


0.230 


4.48 


0.23 


4096 


3.0 


20 


10.0 


13.0 


0.066 


0.074 


0.235 


6.39 


0.17 


4096 


3.0 


50 


9.8 


13.9 


0.075 


0.090 


0.190 


8.84 


0.24 


4096 


3.0 


100 


14.1 


16.5 


0.062 


0.072 


0.274 


34.05 


0.09 
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09 
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77 
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4096 


7.0 


2 


20.8 


51.8 


0.035 


0.045 


0.041 


0.03 


0.52 


4096 


7.0 


2 


17.4 


48.4 


0.036 


0.046 


0.053 


0.06 


0.61 
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o 
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^9 9 
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U. 11U 
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4096 


7.0 


5 


14.0 


53.4 


0.063 


0.078 


0.089 


0.23 


0.69 


4096 


7.0 


10 


3.3 


48.0 


0.100 


0.121 


0.132 


0.61 


0.95 


4096 


7.0 


10 


7.6 


45.8 


0.096 


0.117 


0.126 


0.57 


0.84 


4096 


7.0 


20 


2.0 


41.4 


0.132 


0.158 


0.132 


1.24 


0.98 


4096 


7.0 


20 


6.6 


40.2 


0.097 


0.118 


0.130 


1.44 


0.83 


4096 


7.0 


50 


7.5 


44.0 


0.113 


0.141 


0.163 


4.24 


0.84 


4096 


7.0 


100 


3.4 


55.1 


0.147 


0.189 


0.203 


11.08 


0.93 


4096 


11.0 





0.0 


28.7 


0.019 


0.028 


0.034 


0.02 


1.00 


4096 


11.0 





0.0 


32.2 


0.021 


0.030 


0.031 


0.02 


1.00 


4096 


11.0 


10 


0.0 


25.1 


0.082 


0.103 


0.121 


0.54 


1.00 


4096 


11.0 


10 


0.0 


30.9 


0.079 


0.097 


0.122 


0.56 


1.00 


4096 


11.0 


20 


0.0 


23.8 


0.074 


0.094 


0.128 


1.08 


1.00 


4096 


11.0 


20 


0.0 


26.7 


0.082 


0.102 


0.137 


1.03 


1.00 


4096 


11.0 


50 


0.0 


28.4 


0.080 


0.101 


0.158 


2.86 


1.00 


4096 


11.0 


100 


0.0 


32.4 


0.082 


0.107 


0.179 


11.93 
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16384 


3.0 


10 


9.4 


9.6 


0.034 


0.038 


0.080 


6.05 


0.02 


16384 


3.0 


20 


8.5 


9.2 


0.036 


0.043 


0.241 


6.75 


0.05 


16386 


7.0 


10 


7.2 


32.4 


0.069 


0.086 


0.097 


1.16 


0.76 


16384 


7.0 


20 


6.2 


28.5 


0.076 


0.096 


0.095 


2.53 


0.81 


16384 


11.0 


10 


0.0 


20.1 


0.056 


0.074 


0.102 


0.78 


1.00 


16384 


11.0 


20 


0.0 


19.3 


0.069 


0.090 


0.105 


1.77 


1.00 


512 


7.0 


10 


1.8 


83.9 


0.137 


0.153 


0.196 


0.35 


0.99 (averaged + tidal cut-off) 


1024 


7.0 


10 


3.5 


75.4 


0.124 


0.142 


0.172 


0.39 


0.96 (averaged + tidal cut-off) 


2048 


7.0 


10 


4.7 


67.2 


0.111 


0.132 


0.146 


0.53 


0.94 (averaged + tidal cut-off) 


4096 


7.0 


10 


3.5 


56.3 


0.095 


0.118 


0.129 


0.64 


0.94 (tidal cut-off) 


8192 


7.0 


10 


10.9 


50.4 


0.077 


0.096 


0.108 


0.74 


0.78 (tidal cut-off) 


16384 


7.0 


10 


7.6 


41.2 


0.065 


0.084 


0.088 


0.88 


0.83 (tidal cut-off) 


4096 


3.0 


10 


12.9 


17.0 


0.061 


0.068 


0.275 


2.30 


0.22 (tidal cut-off) 
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3.0 


10 


13.5 
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5.53 


0.09 (tidal cut-off) 


16384 


3.0 


10 


14.2 


15.3 


0.041 


0.046 


0.272 


8.31 


0.07 (tidal cut-off) 


512 


7.0 


10 


1.6 


N/A 


0.141 


0.157 


0.193 


0.32 


1.00 (averaged, isolated) 


1024 


7.0 


10 


3.2 


N/A 


0.131 


0.152 


0.155 


0.40 


0.99 (averaged, isolated) 


2048 


7.0 


10 


4.9 


N/A 


0.112 


0.134 


0.128 


0.50 


0.99 (averaged, isolated) 


4096 


7.0 


10 


5.6 


N/A 


0.104 


0.128 


0.107 


0.59 


0.98 (averaged [5 sim.], isolated) 


8192 


7.0 


10 


9.7 


N/A 


0.082 


0.104 


0.094 


0.64 


0.97 (isolated) 



